score.easy <- function (my.inds, cols = 1, n.inds = NULL, panel = NULL, thresh = NULL,
                        shift = 0.8, ladder, channel.ladder = NULL, ploidy = 2, ci.upp = 1.96,
                        ci.low = 1.96, dev = 50, left.cond = c(0.6, 3), right.cond = 0.35,
                        warn = FALSE, window = 0.5, init.thresh = 200, ladd.init.thresh = 200,
                        method = "cicor", env = parent.frame(), plotting = TRUE, pref = 2, plotfile = "test.png")
{
  if (length(n.inds) > length(my.inds)) {
    print(paste("Hey! you are trying to examine more individuals than the ones you actually read? You selected in 'n.inds' argument",
                length(n.inds), "individuals but you only provided",
                length(my.inds), " individuals. Please select a number of individuals smaller or same size than the ones contained in 'my.inds' argument"))
    stop
  } else {
    #cat(paste("\n1) You have used a shift of", shift,
    #          "base pairs. All peaks at that distance from the tallest peak will be ignored and be considered noise.",
    #          " \n2) In addition the window used is", window, ". Which means that all peaks closer by that distance to panel peaks will be accounted as peaks.",
    #          " \n3) Remember using the get.scores() function to extract the results from this output. \n\n"))
  }
  if (method == "ci") {
    print(paste("Please make sure you have used the same 'dev' value you found convenient for your ladder detection or probably your call will not match"))
  }
  if (is.null(channel.ladder)) {
    channel.ladder <- dim(my.inds[[1]])[2]
  } else {
    channel.ladder <- channel.ladder
  }
  if (dim(my.inds[[1]])[2] < channel.ladder) {
    print(paste("ERROR MY FRIEND!! you have indicated an argument channel.ladder=5, but your data contains less channels/colors"))
    stop
  }
  if (is.null(n.inds)) {
    n.inds <- c(1:length(my.inds))
  } else {
    n.inds <- n.inds
  }
  if (is.null(thresh)) {
    thresh <- rep(list(c(1, 1, 1, 1, 1)), length(my.inds))
  } else {
    thresh <- thresh
  }
  count <- 0
  tot <- length(n.inds)
  #pb <- txtProgressBar(style = 3)
  #setTxtProgressBar(pb, 0)
  my.inds2 <- list(NA)
  thresh2 <- list(NA)
  for (i in 1:length(n.inds)) {
    count <- count + 1
    v1 <- n.inds[i]
    my.inds2[[i]] <- my.inds[[v1]]
    names(my.inds2)[i] <- names(my.inds)[i]
    #setTxtProgressBar(pb, (count/tot) * 0.25)
  }
  ncfp <- c("COLOR 1", "COLOR 2", "COLOR 3", "COLOR 4", "COLOR 5",
            "COLOR 6")
  cfp <- c("cornflowerblue", "chartreuse4", "gold2", "red",
           "orange", "purple")
  col.list <- list(NA)
  att1 <- numeric()
  list.data <- list(NA)
  if (exists("list.data.covarrubias")) {
    list.data <- env$list.data.covarrubias
  } else {
    list.ladders <- lapply(my.inds2, function(x) {
      y <- x[, channel.ladder]
      return(y)
    })
    # import find.ladder
    list.data <- lapply(list.ladders, find.ladder, ladder = ladder,
                        ci.upp = ci.upp, ci.low = ci.low, draw = FALSE, dev = dev,
                        warn = warn, method = method, init.thresh = ladd.init.thresh)
  }
  list.models <- lapply(list.data, function(da) {
    y <- da[[3]]
    x <- da[[1]]
    mod <- lm(y ~ I(x) + I(x^2) + I(x^3) + I(x^4) + I(x^5),
              data = da)
    return(mod)
  })
  list.models.inv <- lapply(list.data, function(da) {
    x <- da[[3]]
    y <- da[[1]]
    mod <- lm(y ~ x, data = da)
    return(mod)
  })
  xx <- lapply(my.inds2, function(x, cols) {
    1:length(x[, cols])
  }, cols = cols)
  newxx <- numeric()
  newyy <- numeric()
  new.whole.data <- list(NA)
  for (h in 1:length(xx)) {
    h1 <- n.inds[h]
    count <- count + 1
    newxx <- as.vector(predict(list.models[[h1]], newdata = data.frame(x = xx[[h]])))
    newyy <- my.inds2[[h]][, cols]
    new.whole.data[[h]] <- list(xx = newxx, yy = newyy)
    #setTxtProgressBar(pb, (count/tot) * 0.25)
  }
  top <- max(unlist(lapply(new.whole.data, function(x) {
    max(x$yy)
  })))
  bott <- min(unlist(lapply(new.whole.data, function(x) {
    min(x$yy)
  })))
  list.weis <- list(NA)
  lower.bounds <- numeric()
  for (k in 1:length(my.inds2)) {
    newtt <- init.thresh
    lower.bounds[k] <- newtt
    plant <- big.peaks.col(new.whole.data[[k]]$yy, newtt)
    plant$wei <- new.whole.data[[k]]$xx[plant$pos]
    # import separate
    plant <- separate(plant, shift, type = "bp")
    list.weis[[k]] <- plant
    if (plotting == TRUE) {
      count <- count + 1
    } else {
      count <- count + 2
    }
    #setTxtProgressBar(pb, (count/tot) * 0.25)
  }
  list.weis <- lapply(list.weis, function(x) {
    x$wei <- round(x$wei, digits = 4)
    return(x)
  })
  names(list.weis) <- names(my.inds2)
  if (length(panel) > 0) {
    # import reals
    list.weis <- lapply(list.weis, reals, panel = panel,
                        shi = shift, ploidy = ploidy, left.cond = left.cond,
                        right.cond = right.cond, window = window)
    # import homo.panel
    list.weis2 <- lapply(list.weis, FUN = homo.panel, panel = panel,
                         window = window)
  } else {
    list.weis2 <- list.weis
  }
  if (plotting == TRUE) {
    png(plotfile, width = 2000, height = 800)
    par(mfrow = c(pref,1))
    #layout(matrix(1:pref, pref, 1))
    marker.list <- list()
    if (cols == 1) {
      marker.list <- B.marker
    } else if (cols == 2) {
      marker.list <- G.marker
    } else if (cols == 3) {
      marker.list <- Y.marker
    } else if (cols == 4) {
      marker.list <- R.marker
    }
    #str(marker.list)
    if (length(panel) > 0) {
      xm <- round(min(panel, na.rm = TRUE) - 10, digits = 0)
      xl <- round(max(panel, na.rm = TRUE) + 10, digits = 0)
    } else {
      xm <- 0
      xl <- max(ladder)
    }
    for (g in 1:length(n.inds)) {
      hh4 <- n.inds[g]
      if (length(which(new.whole.data[[g]]$xx > xm & new.whole.data[[g]]$xx < xl)) > 0) {
        mylim <- max(
            new.whole.data[[g]]$yy[which(new.whole.data[[g]]$xx > xm &
            new.whole.data[[g]]$xx < xl)], na.rm = TRUE)
        + 100
      } else {
        mylim = 1000
      }
      if (is.infinite(mylim)) {
        mylim = 1000
      }
      plot(new.whole.data[[g]]$xx, new.whole.data[[g]]$yy,
           type = "l", col = cfp[cols], xaxt = "n", xlim = c(xm, xl),
           ylim = c(-200, mylim+1000), ylab = "",
           #main = paste(ncfp[cols], "plant", hh4),
           xlab = names(list.models)[hh4],
           lwd = 2, las = 2)
      axis(1, at = c(xm:xl), labels = xm:xl, cex.axis = 0.8)
      #rect(xleft = (list.weis2[[g]]$wei - window),
      #     ybottom = (bott - 200), xright = (list.weis2[[g]]$wei + window),
      #    ytop = (top + 1000), col = transp("lightpink", 0.3), border = NA)
      #str(list.weis2[[g]]$wei - window)
      for (marker in names(marker.list)) {
        xleft = marker.list[[marker]][1]
        xright = marker.list[[marker]][[2]]
        ybottom = mylim+700
        ytop = mylim+1000
        rect(xleft = xleft, ybottom = ybottom, xright = xright, ytop = ytop)
        center <- c(mean(c(xleft,xright)), mean(c(ybottom, ytop)))
        #str(center)
        text(center[1],center[2],labels = marker)
        #text(xright,ytop,labels=marker)
      }
      #abline(v = list.weis[[g]]$wei, lty = 3, col = "blue", cex = 0.5)
      #print(list.weis[[g]]$wei)
      #abline(v = list.weis2[[g]]$wei, lty = 3, col = "red", cex = 0.5)
      #abline(h = lower.bounds[g], lty = 2, col = "chocolate", cex = 0.5)
      #legend("topright", legend = c("Peak found", "Panel peak", "Panel window", "Minimum Detected"),
      #       col = c("blue", "red", transp("lightpink", 0.3), "chocolate"),
      #       bty = "n", lty = c(3, 3, 1, 3), lwd = c(1, 1, 3, 1), cex = 0.75)
      count <- count + 1
      #setTxtProgressBar(pb, (count/tot) * 0.25)
    }
    dev.off()
  }
  #close(pb)
  #return(list.weis2)
  #return(list(newxx,newyy))
  #return(new.whole.data)
  #return(top)
  #return(1)
}
